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ABSTRACT 

New predictive closure models for turbulent free shear flows are presented in this paper. They 
are based on an instability wave description of the dominant large-scale structures in these flows 
using a quasi-linear theory. Three models have been developed to study the structural dynamics 
of turbulent motions of different scales in free shear flows. The local characteristics of the large- 
scale motions are described using linear theory. Their amplitude is determined from an energy 
integral analysis. The models have been applied to the study of an incompressible free mixing 
layer. In all cases, predictions are made for the development of the mean flow field. In the last 
model , predictions of the time-dependent motion of the large-scale structure of the mixing region 
are made. The predictions show good agreement with experimental observations. 

INTRODUCTION 

Though the presence and importance of large-scale coherent structures to the mixing process 
in free shear flows has been recognized for many years, turbulence models thai incorporate these 
observations have been very limited. The use of direct numerical or large eddy simulations provide 
a detailed prediction of the large-scale motions in low and high Reynolds number turbulent flows 
respectively. But these predictions are computationally expensive and are still limited in general 
to simple boundary conditions. The present model makes use of experimental observations in 
excited turbulent flows or conditional sampling measurements to provide a simple model of the 
large scale motions which is computationally inexpensive. 

Most current turbulent flow calculations for practical applications use the long time-averaged 
Navier-Stokes equations. Turbulence models are needed to evaluate the unknown correlation 
terms, the Reynolds stresses, that appear when the statistical averaging process is applied to 
the nonlinear convective terms in the equations. This is the closure problem. There are closure 
models of various orders that have been proposed. These models are usually based on the notion 
that the high-order moments of fluctuations can be represented reasonably well as functionals of 
moments of lower order. Work in this regard is voluminous and will not be elaborated on here. 
Some models are quite successful and have become very popular in engineering flow calculations. 
However, they usually involve a large number of model constants determined by comparison with 
experimental data. Thus, these models are not entirely predictive but, in some ways, represent a 
sophisticated correlation of experimental da^a. 
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The present models are based on observations of large-scale coherent structures in free mixing 
layers. Brown and Roshko (ref.l), among others, observed that these orderly motions dominate 
the dynamics and the structure of free shear flows like wakes, jets and mixing layers. The 
structures appear in both low- and high-speed flows. They have also been observed in many flow 
geometries. 

This paper is concerned with new, predictive turbulence models for free shear flows. The models 
simulate the propagating large-scale structures as spatially travelling instability waves. In this 
paper, we focus on the validation of the wave models as well as a determination of their limitations. 
Predictions are made for a two-dimensional incompressible free mixing layer. This will provide 
guidelines for applications of the models to more complex configurations. 

THE WAVE MODELS 

The wave models presented here are used to make a direct calculation of the large-scale, 
characteristic structures. The fundamental idea is that the large-scale structures may be mod- 
eled using a quasi-linear theory. The local characteristics of these structures may be described by 
linear instability theory. This has been demonstrated by the experiments of Gaster, Kit and Wyg- 
nanski (ref.2) and Petersen and Samet (ref.3), among others. In their experiments they compared 
predictions of the amplitude and phase of the axial velocity fluctuations, based on linear stability 
theory, with phase-averaged measurements in an excited shear layer and a jet. The agreement 
between predictions and experiment was very good though only normalized distributions of am- 
plitude and phase, not the absolute amplitude, were predicted. This close agreement between 
the predictions of linear stability theory and the properties of the large-scale coherent sturctures 
has formed the basis for theories of turbulent mixing and supersonic jet noise generation and 
radiation. For example, Tam and Morris (ref.4) and Tam and Burton (ref.5, ref.6) predicted the 
noise radiation from instability waves in supersonic shear layers and jets and obtained very good 
agreement with experiment. Their analyses showed that the behavior of the large-scale distur- 
bances could be modeled satisfactorily using a quasi-linear theory, even though the waves were 
not infinitesimal in magnitude. However, an important element of these calculations, the velocity 
profiles of the mean flow, that are needed for the linear stability calculations, are obtained from 
experiments. Their approaches provide a closure, but are not predictive. The models proposed 
here establish a complete closure model using a simple quasi-linear theory for the large-scale 
motion. In the present model both the mean flow and the time-dependent turbulent motions at 
the large-scale are obtained simultaneously as a solution. 

Analysis 

In the present analysis, the turbulent motion is decomposed into three parts, 

h = Fi + U + /'. ( 1 ) 

The fluctuation with respect to the mean flow,Iq, is separated into a component representing the 
dominant large-scale motion, fi , and one representing small-scale fluctuations, ff The mean flow 
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component is obtained by long time-averaging its instantaneous value 


fi = Fi = — / fid t. 

1 1 Jo 

A short time-average is defined by 


< fi > — F{ + fi 



( 2 ) 

(3) 


where T 2 is much smaller than T\, but much larger than the characteristic time scale of the 
background small-scale fluctuation. 


The governing equations for the mean flow can be obtained by substituting flow properties of the 
form of equation (1) into the Navier-Stokes equations and long time-averaging the equations. 
The governing equations for the mean flow are 
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where the interactions between motions of disparate scales are assumed to be negligible. The 
boundary-layer approximation is used to simplify further the governing equations in the present 
case. The resulting equations are 
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where the small-scale normal stress terms have been neglected. On the basis of the experimental 
observations the local characteristics of the large-scale structures are described by the equation 
of inviscid hydrodynamic stability. The equations of motion for the large-scale fluctuations are 
linearized. Solutions are sought in the form: 


{u,v,p} = A(x) dt{[(f>(y),ip(y),<p(y)} exp[i{ax- wt)]}. (8) 

where the amplitude appears as a parameter in the local problem and is determined separately. 
This weakly non-linear approach is usually referred to as a “wave envelope” method. The resulting 
equation for the fluctuating velocity in the cross-stream direction, y, is the Rayleigh equation. 

( aU — u) ) {if" - a 2 if) - aU"if = 0. (9) 

where the main stream is in the x-direction of a Cartesian co-ordinate system and ( )" denotes 
d 2 /(dy) 2 . In deriving this equation, it is assumed that the mean flow is locally parallel in the 
direction of the main stream. 
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The amplitude function, A(x), can be determined from the kinetic energy equation for the large- 
scale motions, 
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where k = 1/2 u t u r . At the large scale, the viscous terms are negligible. The production 
terms on the right hand side are responsible for transferring energy from the mean flow to the 
coherent turbulent fluctuations. Energy is subsequently extracted from the large-scale motion 
and dissipated at the high frequency end of the wave number spectrum. The terms containing the 
residual stress tensor, — < u'u'- >, describe the draining of energy from the waves. These terms 
are of crucial importance in determining the wave amplitude. Little experimental information, 
however, is available regarding these stresses. In the present paper several approaches have been 
taken to model this energy transfer mechanism. Initially we assume that the rate of energy 
dissipation is proportional to 


where l and u are the characteristic length and velocity scales of the large-scale motions. This 
model assumes that the turbulence is in an equilibrium state for the small-scale fluctuations, in 
which the rate at which energy is transferred from the large scales is equal to the rate at which 
energy is dissipated. The net effect of these terms may thus be modeled by 


r k 3 / 2 

i — r~ 


( 12 ) 


where C\ is a model constant. An equation for the amplitude function may be obtained by 
substituting the wave form expressions, equation (8), into the wave kinetic energy equation, 
equation (10), and integrating it with respect to y. The resulting equation is 
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where an asterisk denotes the complex conjugate. 
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A complete simulation of the large-scale turbulence spectrum would require the inclusion of a 
broad range of frequency and spanwise wave number components. This was accomplished in the 
local solution of Tam and Chen (ref. 7) and the the integral model of Morris and Giridharan (ref.8). 
However, it can be shown that for a wide range of frequencies around the least stable mode the 
Reynolds stress distribution does not vary appreciably. Since sensible and manageable models of 
maximum efficiency are sought, instead of including all the unstable waves, it is assumed here 
that the waves of the least stable modes are most effective in extracting energy from the mean 
flow and are used to describe the overall properties of the coherent structures. 

The contribution of the small-scale Reynolds stress gradients in equation (7) has yet to be deter- 
mined. In the first model described below this contribution is given by a simple eddy viscosity 
model. Thus, this model, referred to as Model I, requiries the specification of two model con- 
stants: one to determine the rate of energy transfer from large to small scales, equation (12), and 
one for the eddy viscosity model. However, the turbulence models of Tam and Chen (ref. 7) and 
of Morris and Giridharan (ref.8) suggest that the small scales need not play a direct role in the 
development of the mean flow. Thus in Model II we neglect the contributions of the small-scale 
Reynolds stresses. Both of these models predict the average behavior of the shear layer. 

A third model, Model III, will also be examined. This model simulates the time-dependent 
motion, at the large scale, associated with the passage of a train of large-scale structures. Ex- 
perimental observations suggested that, even if initially there exists a continuous spectrum of 
infinitesimal disturbances upstream of the splitter plate, a disturbance emerges dominating over 
other neighboring perturbations in the early stages of the flow development. As the flow evolves, 
however, there are continuous shifts of the dominant component toward lower and lower fre- 
quencies. In fact, the growth of an initially small periodic disturbance is often followed by the 
development of subharmonics. In numerical simulations, however, the initial conditions can be 
conceived in a much simplier way. Instead of monitoring the disturbances in the complete initial 
continuous spectrum, a hierarchy of disturbances made up of the initially most unstable mode, 
according to linear theory, and its subharmonics may be considered. This reflects the “subhar- 
monic evolution model” proposed by Ho and Huang(ref.9). The unsteady turbulent fluctuations 
of large scale are thus described by the superposition of the instability waves in this hierarchy. 
This enables the time-dependent flow field at the large scale to be simulated. 

The mean flow and the local shapes of the large-scale structures or instability waves are governed 
by, as stated earlier, the thin-shear-layer and Rayleigh equations, respectively. The solution 
methods for these equations are the same as in Models I and II. The equation for the amplitude 
function, however, has to be modified. Firstly, it is assumed that interactions between harmonics 
are negligible, as there is sufficient phase jitter in the unexcited shear layer. In addition, the details 
of the process of energy transfer from large to small scales is not modeled explicitly. At the axial 
location where a given instability wave saturates, or begins to decay, the energy is immediately 
removed from the system. Thus there is no need to specify either a constant associated with 
the energy transfer process or the effects of the interaction between the small-scales motion and 
the mean flow. There are no empirical constants. The amplitude of each instability wave during 
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the unstable or growing region is determined from equation (13) with the energy transfer terms 
neglected. 

A visualization of the unsteady flow field predicted by Model III is made by means of instantaneous 
flow velocity vector plots and streaklines. The streaklines are produced by injecting passive marker 
particles at the initial location, x = xo at various points across the shear layer. The positions of 
these particles at subsequent times can be calculated using the equations: 

= U [ x{t) , y{t) } + u [ x(t) , y(t ) , t ] (17) 

at 

-^y(t) = V [ x(t) , y(t) ) ] + v l x{t) , y(t) , t ], (18) 

x(0) = x 0 , 2/fc (0) =y 0 , k = 1 
Particles thus move at each time step according to the local instantaneous velocity field. 

Numerical Procedure 


and 


with 


The boundary-layer approximation renders the system of equations governing the mean flow 
parabolic. A fourth order Keller-Box scheme is applied to solve the resulting equations. The 
equation for the instability wave, which is the Rayleigh equation in the present formulation, has 
been solved using various methods, including a traditional shooting, two spectral and a finite 
difference methods. For spatial instability, the system of equations generated by the global 
approximations of the Rayleigh equation forms an eigenvalue problem which is nonlinear in its 
parameter, the wave number. It may be solved using the Linear Companion Matrix method or a 
method based on matrix factorization, Bridges and Morris (ref.10). Details of the various solution 
schemes can be found in Liou and Morris (refill). The Rayleigh equation and the equations for 
the mean flow are solved iteratively at each streamwise location. The convergence criterion for 
the iterations is 
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where ej is a small number and M is the total number of grid points at each streamwise location. 
The amplitudes of the waves are calculated explicitly using a fourth order Runge-Kutta method 
applied to the wave energy equation (13). 


RESULTS AND DISCUSSIONS 


The models have been tested in an incompressible free mixing layer. The flow is sketched in 
figure (1). A hyperbolic tangent distribution is taken as the shape of the initial streamwise mean 
velocity, U(xo,y), i.e., 


U ( x 0 , y ) 



tanh ( 30 y )). 


( 20 ) 
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Figure 1. Sketch of a free mixing layer. 
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Figure 3. -ut> in the self-similar region, fc , 
Chebyshev Collocation (N=ll); A , Cheby- 
shev Collocation (N=19); o , finite difference 
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Figure 5. Growth of the mixing layer using 
Model I. 



Figure 2. Physical domain. 
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The initial cross-stream mean velocity, V (xq , y), can be set to a small number or zero. The 
boundary conditions for the mean flow are 

U ( x , y u (x) ) = 1.0, U ( x , yi(x) ) = 0.0, V ( x , y u (x) ) = 0.0. (21) 

where y u {x) and yi{x) are the upper and the lower boundaries of the physical domain shown in 
figure (2). As a test of the ability of the instability wave model to describe large-scale structures 
and the associated Reynolds stresses, the model was first applied in the self-similar region of 
the flow with a mean velocity profile from Patel (ref. 12). Figure (3) shows the calculated and 
experimental Reynolds shear stress distributions. Calculated results using a traditional shooting 
method compare favorably with that using global approximations of various order. Note that all 
the calculated results have been normalized by the peak experimental value. The discrepancy at 
the low-speed side of the layer suggested that the momentum exchanges due to the small-scale 
turbulent motions might not be negligible in this region. It should be noted that this negative 
value of Reynolds shear stress disappears for small values of velocity in the lower stream. The 
structures obviously contribute negative shear stress at the low-speed edge of the flow. This 
counter-gradient transport of momentum gives negative energy production in this region. A 
similar phenomenon was observed experimentally by Komori and Ueda (ref. 13) in the self-similar 
region of a jet. In fact, regions of negative shear stress can be easily observed if the large-scale 
structures are excited artificially, for example, see Wygnanski, Oster and Fiedler (ref.14). This 
counter-gradient momentum transfer decelerates and subsequently reverses the flow on the low- 
speed. side of the mixing layer as the shear layer develops. As noted above, Model I proposes that 
a contribution from the small-scale Reynolds stresses is required to describe the total turbulent 
forces that determine the development of the mean flow. Thus, in Model I, a simple eddy-viscosity 
model is used to model the small-scale Reynolds shear stress, 


— u'v' 
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(22) 


The model introduces a new parameter, C 2 . Latigo (ref.15) argued that the turbulent shear stress 
contributed by the small-scale motions is about a half of the total shear stress. An estimate of 
C 2 based on the value that is used in the classical eddy-viscosity models is then obtained. In 
addition, the force terms associated with the large-scale normal stresses in the mean momentum 
equations are also retained, since they are found to be of the same order as the other Reynolds 
stress gradient terms on the low-speed edge of the shear layer. The normal stresses associated 
with the large-scale structures can be calculated directly by the wave models and involve no 
further empiricism. 


In the numerical calculations, the local solution of the Rayleigh equation is found to be time- 
consuming. To accelerate the axial marching an adaptive grid has been devised. The grid size in 
the cross-stream direction in the transformed domain are fixed. The axial step sizes are chosen 
such that the convergence indices of the first iteration at a downstream station are greater than 
a fixed number £ 2 - 
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The grids are found to cluster in a region where there are large changes of flow properties, for 
example, when the flow is developing initially. 

The initial wave amplitude represents the initial strength of the instability waves or large-scale 
motions for which there are no quantitative experimental measurements. From a sequence of 
numerical experiments, however, it is found that flows with relatively strong initial amplitude 
saturate early. Subsequently the flow develops in a similar manner; only the virtual origin of 
mixing is changed. The initial amplitudes for the cases presented in this paper are fixed at 10~ 4 . 
The model constant Cb of the energy transfer term in the wave kinetic energy equation is taken 
from a conventional Prandtl energy model, Launder, et. al.(ref.l6). It is found that its value has 
no significant influence on the results of the mixing layer calculations. 

First we consider Model I in which turbulent forces associated with the wave shear and normal 
stresses as well as the small-scale motions are considered. Figure (4) shows the axial forces acting 
on fluid elements across the layer at various axial stations. The calculated rate of growth of the 
layer agrees well with the value that is an average taken over various experiments. The averaged 
experimental value is denoted by the straight line in figure (5). 8 is the distance between the 
points where the local mean velocity is 0.9 and 0.1 of the main stream mean velocity, U. Figure 
(6) shows the predicted mean velocities at a sequence of downstream stations. The agreement 
between the predictions and experimental data is good except near the low-speed side of the 
shear layer. Note that for a free mixing layer, the accuracy of the measured mean flow data in 
this region is poor due to the rapid variations in the instantaneous flow direction. As can be 
seen from figure (7), which shows the shear stress distributions across the mixing layer at various 
stations, the sum of the shear stresses from the large-scale and the small-scale motions agrees well 
with experimental data. The experimental measurements are the long time-averaged correlations 
of the turbulent fluctuations. The amplitude of the large-scale fluctuations is plotted in figure 
(8). The large-scale structures extract energy from the mean flow and strengthen as the flow 
develops. However, energy is also being transferred to the small scales. The final equilibrium of 
the large-scale motions amplitude is reached when the energy gained from the mean flow balances 
the energy lost to the small scales. 

In Model I, the small-scale motions play a direct role in the momentum transport process. How- 
ever, it could be argued that the large-scale structures dominate dynamically the development 
of free shear flows. Therefore, in Model II only the fluctuations at the large scale are included. 
This eliminates the need to specify a model constant in the eddy-viscosity model, equation (22). 
It can be seen from figure (9) that the forces associated with the large-scale normal stresses are 
apparently able to counter-balance the decelerating effects of the wave shear stress gradients. 
The predicted mean velocity profiles are presented in figure (10). It shows that the mean flow 
can be satisfactorily predicted by modeling only the dominant large-scale structures. However, 
as shown in figure (11), the predicted shear stress distributions do not match the total shear 
stress distributions measured by Patel (ref. 12). This difference does not mean necessarily that 
the smll-scale stresses should be included. It must be remembered that the present model simu- 
lates the entire large-scale spectrum with a single frequency wave that is locally most unstable. 
Tam and Chen (ref. 7), in their local model, included a broad range of instability waves and found 
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Figure 11. Distributions of wave shear stress 
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Figure 8. Evolution of the large-scale struc- 
ture amplitude using Model I. 
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Figure 10. Mean velocity profiles using Model 

II at , £=0.63; , 2.96; , 7.36; 

, 9.36; , 11.36. a , Patel. 



V 
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good agreement with experiment, without the inclusion of contributions from the small scales. 
This point requires further investigation. The evolution of the large-scale amplitude using Model 
II follows a similar behavior to that shown in figure (8). Once again an equilibrium condition 
is reached where the rate of energy transfer from the mean flow to the large-scale structures 
balances the rate at which energy is lost by the structures to small scales for eventual dissipation. 

Model III simulates the time-dependent motion, at the large scale, associated with the passage of 
a train of large-scale structures. As the flow develops axially, these hydrodynamic waves become 
damped because of the growth of the shear layer. Since it is assumed that energy associated 
with a given wave is removed immediately it becomes neutral, there is no need to obtain damped 
inviscid solutions by analytic continuation in the complex plane, Tam and Morris (ref. 17). 

In the preliminary calculations, it was found that an abnormality in the mean velocity distribu- 
tions appeared near the critical points of saturating waves. Also, most of the shear layer growth 
occurred on the low speed side of the layer. This gives a non-monotonic velocity distribution 
near the critical layers of saturating waves and another inflection point appears. Saturating waves 
thus have to be removed before they become neutral during the axial marching. Wygnanski and 
Petersen (ref. 18) suggested that this abnormality is due to nonlinearities. Composite expansion 
techniques have been applied to investigate the effect of critical-layer nonlinearity, for example, 
see Goldstein and Leib (ref. 19) and Goldstein and Hultgren (ref.20). Another approach to resolve 
this issue is to include viscous effects; that is to solve the Orr-Sommerfeld equation. Since the 
present investigation is directed toward developing simple turbulence models, instead of includ- 
ing other computationally expensive approaches, the effects of the critical point is accounted for 
by incorporating a small amount of eddy-viscosity in the analysis of the mean flow. With this 
modification the mean velocity distributions predicted at several downstream stations are shown 
in figure (12). They are compared with Patel’s data using a similarity scale. The additional 
mixing at the fine scale is diffusive and able to smooth out humps in the flow. In the present 
calculation, the extra mixedness provided accounts for about 10% of the amount of turbulent 
momentum exchange that is suggested by conventional models. There are six waves in the hier- 
archy in this calculation. Since waves are removed successively during the axial marching, the 
number of waves included depends on the distance the calculation is to be carried downstream. 
There are some differences between the calculated results and Patel’s measurements in figure (12). 
However, similar characteristics of mean velocity distributions have been reported by Wygnanski 
et. at. (ref. 14), among others, in which free mixing layers are excited externally. Figure (13) 
shows the development of the wave amplitudes. The additional small-scale mixing increases the 
initial growth of the layer so that the fundamental mode is removed at a lower amplitude than 
its subharmonics before its amplitude reaches equilibrium level. The axial width of the layer is 
shown in figure (14) and is compared with the prediction using Model I. The presence of this 
stepwise evolution is characteristic of excited flows and would be smoothed out if many waves 
with slightly different amplitudes and frequencies were included. 

Figure (15) shows the transient solutions of the velocity field of the turbulent free mixing layer 
depicted by Model III in a reference of frame moving at the phase speed of the fundamental wave. 
Dominant large vorticies can be clearly seen. The downstream development of the neighboring 
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Figure 15. Velocity vector plots at (a) t=1.5, 
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vortical structures can also be observed. This can be assisted by streakline plots shown in figure 
(16). The roll-up of vorticies into larger vortex-like structures can be observed clearly. The initial 
roll-up is dominated by the fundamental mode. As time progresses, the initial structures convect 
downstream and roll around each other. These regions of concentrated vorticity then form a single 
large structure. As the passive particles travel downstream, their motion becomes dominated by 
lower subharmonics. Vortex-like structures of increasing scale are formed. Subsequently, the 
rolling process between two adjacent structures repeats as the flow develops further downstream. 
Careful examination of the figures shows also how the structures are convecting downstream as 
they form and roll. Large tongues of unmixed fluid are swept across the layer and reach the 
opposite side of the layer as observed by Brown and Roshko (ref.l). The engulfed fluid elements 
from the two sides of the layer mix and are drawn into the leading and trailing vortices when 
passing through the high-strain braid region between the vortices. This provides the environment 
for further fine-scale mixing. 


SUMMARY 

Three models based on a quasi-linear theory, that describes the dynamics of the dominant 
large-scale structures in a free mixing layer have been presented. The closure schemes incorpo- 
rating the models are able predict the development of the turbulent free mixing layer accurately, 
even though they contain some assumptions and simplifications. The transient turbulent motions 
at the large scale in the layer mapped out using Model III possess many features that are apparent 
in flow visualization experiments, such as the convective nature of the large-scale structures, the 
large-scale transport of unmixed fluid elements and the roll-up of vortices. The models involve 
less empiricism than most conventional models. Since large-scale coherent structures appear 
also in shear flows of other geometries, the closure schemes presented here should be applicable 
to those cases as well. It is hoped that these models, which originate from observed physical 
phenomena, will provide efficient tools to model other free shear flows. 
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